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Abstract 



The processes influencing animal movement and resource selection are 
complex and varied. Past efforts to model changing behavior over time used 
Bayesian statistical models with variable parameter space, such as reversible- 
jump Markov chain Monte Carlo approaches, which are computationally 
demanding and inaccessible to many practitioners. We present a continuous- 
time discrete-space (CTDS) model of animal movement that can be fit using 
standard generalized linear modeling (GLM) methods. This CTDS approach 
allows for the joint modeling of location-based as well as directional drivers 
of movement. Changing behavior over time is modeled using a varying- 
coefficient framework which maintains the computational simplicity of a 
GLM approach, and variable selection is accomplished using a group lasso 
penalty. We apply our approach to a study of two mountain lions (Puma 
concolor) in Colorado, USA. 

Keywords: Agent-Based Model; Animal Movement; Mountain Lion; Mul- 
tiple Imputation; Varying- Coefficient Model. 
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1 Introduction 



Animal telemetry data have been used extensively in recent years to study 
animal movement, space use, and resource selection (e.g., Johnson et al., 
2011; Hanks et al., 2011; Fieberg et al., 2010). The ease with which teleme- 
try data are being collected is increasing, leading to vast increases in the 
number of animals being monitored, as well as the temporal resolution at 
which telemetry locations are obtained (Cagnacci et al., 2010). This com- 
bination can result in huge amounts of telemetry data on a single animal 
population under study. Additionally, the processes driving animal move- 
ment are complex, varied, and changing over time. For example, animal be- 
havior could be driven by the local environment (e.g., Hooten et al., 2010), 
by conspecifics or predator/prey interactions (e.g., Merrill et al., 2010), by 
internal states and needs (e.g., Nathan et al., 2008), or by memory (e.g., 
Van Moorter et al., 2009). The animal's response to each of these drivers of 
movement is also likely to change drastically over time (e.g., Nathan et al., 
2008; Hanks et al., 2011; McClintock et al., 2012). 

The large amount of telemetry data available, even for one animal, and 
the complex behavior displayed in animal movement results in a challenging 
situation for statistical modeling. There is no shortage of existing statistical 
models of animal movement; however, most of these models are compu- 
tationally demanding, and most are inaccessible to the practitioner. For 
example, consider the agent-based model of animal movement of Hooten 
et al. (2010). The agent-based framework is highly flexible, allowing for 
location-based and directional drivers of movement, but is computationally 
expensive. Analyzing the movement path of one animal with hourly loca- 
tions over the course of a week using the approach of Hooten et al. (2010) can 
require computational time on the order of days using standard computing 
resources. The velocity-based framework for modeling animal movement of 
Hanks et al. (2011) allows for time-varying behavior through a changepoint 
model of response to drivers of movement, and is more computationally ef- 
ficient than the approach of Hooten et al. (2010), requiring computational 
time on the order of hours for a similar problem. Similarly, the mechanistic 
state-switching approach of McClintock et al. (2012) allows for time- varying 
behavior through a state-switching approach. These three approaches use 
Bayesian statistical models, and both Hanks et al. (2011) and McClintock 
et al. (2012) allow for time- varying behavior by letting the model parameter 
space vary, either through a reversible- jump Markov chain Monte Carlo ap- 
proach (Green, 1995) or the related birth-death Markov chain Monte Carlo 
approach (Stephens, 2000). These methods can be quite computationally 
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demanding, require the user to tune the algorithm to ensure convergence, 
and can be inaccessible to many practitioners. 

The agent-based model of Hooten et al. (2010) assumes a representation 
of the animal's movement path that is discrete in both space (grid cells) 
and time (fixed time intervals). The velocity-based movement model of 
Hanks et al. (2011) assumes a representation of the movement path that 
is continuous in space and discrete in time. The state-switching model of 
McClintock et al. (2012) assumes a representation of the movement path 
that is discrete in time and continuous in space. 

In this paper, we present a continuous-time, discrete-space (CTDS) model 
for animal movement which allows for flexible modeling of an animal's re- 
sponse to drivers of movement in a computationally efficient framework. 
Instead of a Bayesian approach, we adopt a likelihood-based approach for 
inference, and instead of a state-switching or change-point model for chang- 
ing behavior over time, we adopt a time-varying coefficient model. We also 
allow for variable selection using a lasso penalty. This CTDS approach is 
highly computationally efficient, requiring only minutes or seconds to ana- 
lyze movement paths that would require hours using the approach of Hanks 
et al. (2011) or days using the approach of Hooten et al. (2010), allowing 
the analysis of longer movement paths and more complex behavior than has 
been previously possible. To make this CTDS approach for modeling animal 
movement and resource selection accessible to practitioners, code to imple- 
ment this approach is available online (www.stat.colostate.edu/~hanks) in 
the form of a package for the R statistical computing environment (R De- 
velopment Core Team, 2012) with worked examples. 

In Section 2 Preliminaries, we describe the continuous-time continuous- 
space model of Johnson et al. (2008) which is used to make inference on 
the posterior predictive distribution of an animal's continuous movement 
path, conditioned on observed telemetry locations. We then describe the 
method of multiple imputation (Rubin, 1987) which we use to integrate 
over the uncertainty in the animal's continuous movement path. In Sec- 
tion 3 Continuous-Time Discrete-Space Movement Model, we describe the 
CTDS model for animal movement, and show how inference can be made 
on parameters in this model using standard software for generalized linear 
models (GLMs). In Section 4 Time- Varying Behavior and Variable Selec- 
tion we use a varying-coefficient approach to model changing behavior over 
time, and use a lasso penalty for variable selection. In Section 5 Drivers of 
Animal Movement we discuss modeling potential covariates in the CTDS 
framework. In Section 6 Example: Mountain Lions in Colorado we illus- 
trate our approach through an analysis of mountain lion (Puma concolor) 
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movement in Colorado, USA. Finally, in Section 7 Discussion we discuss 
possible extensions to the CTDS approach. 



2 Preliminaries 

2.1 Continuous-Time Continuous-Space Movement Model 

To model animal movement, we make use of the continuous time correlated 
random walk (CTCRW) model of Johnson et al. (2008) to characterize a dis- 
tribution for the continuous path conditioned on observed telemetry data. 
Let S = {s(t),t = to,t±, . . . , tx} be a collection of time-referenced telemetry 
locations for an animal. If the animal's location and velocity at an arbi- 
trary time t are s(i) and v(i), respectively, then the CTCRW model can be 
specified as follows, ignoring the multivariate notation for simplicity: 



where -0 = [ipi, ip2-, ^3] control the movement and uj(t) is standard Brown- 
ian motion. This model can be discretized and formulated as a state-space 
model, which allows for efficient computation of discretized paths S at arbi- 
trarily fine time intervals via the Kalman filter (Johnson et al., 2008). If a 
Bayesian framework is used for inference on if>, Johnson et al. (2008) shows 
how the posterior predictive distribution of the animal's continuous path 
S can be approximated using importance sampling. We will refer to the 
posterior predictive path distribution as [S|S], where the bracket notation 
'[•]' denotes a probability distribution. 

2.2 Multiple Imputation 

Our general strategy is to construct a model conditioned on the continuous 
path S, and then integrate over the uncertainty in the posterior predictive 
distribution [S|S] (e.g., Hooten et al., 2010; Hanks et al., 2011). If we treat 
the unobserved continuous path S as missing data, then we can make in- 
ference on model parameters using multiple imputation (Rubin, 1987). We 
motivate multiple imputation as posterior predictive inference on the impu- 
tation distribution within a Bayesian framework. Our treatment is similar 
to that of Rubin (1987) and Rubin (1996). 




(1) 



s(t) = s(0) + / v(u)du , 
Jo 
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If we desire posterior predictive inference [6\S] concerning environmen- 
tally relevant movement parameters 0, conditioned on the telemetry data S 
and the posterior predictive path distribution [S|S], then we can write: 

[e\s} = Jje\s][s\s]ds. (2) 

In the multiple imputation literature, the posterior predictive path distri- 
bution [S|S] is called the imputation distribution and specifies a statistical 
model for the missing data S conditioned on the observed data S. We will 
use the CTCRW model of Johnson et al. (2008) as the imputation distri- 
bution [S | S] in our CTDS approach to modeling animal movement. The 
CTCRW model is a mechanistic model of animal movement that has been 
successfully applied to studies of aquatic (Johnson et al., 2008) and ter- 
restrial (Hooten et al., 2010) animals, and can represent a wide range of 
behavior. 

Hooten et al. (2010) and Hanks et al. (2011) use composition sampling 
to obtain samples from the posterior predictive distribution [0\S] in (2) 
by sampling iteratively from [0\S] and [S|S]. Alternately, under the mul- 
tiple imputation framework the posterior distribution [0\S] is assumed to 
be asymptotically Gaussian. The posterior can then be approximated us- 
ing only the posterior predictive mean and variance, which can be obtained 
using conditional mean and variance formulae: 



E(0\S) « f f [6\S][S\S]dSd0 
Je Js 

= j (J^ 0[0\$\d0^j [S|S]dS 

= £7 §|s (E(0\S)) (3) 

and likewise: 

Var(0|S) « £7 §|s (Var(0|S)) + Var g|s (e(0\S)) . (4) 

As the posterior distribution [0\S] converges asymptotically to the sam- 
pling distribution of the maximum likelihood estimate (MLE) of 0, we can 
approximate [0\S] by obtaining the asymptotic sampling distribution of the 
MLE. This allows us to use standard maximum likelihood approaches for 
inference, which can be much more computationally efficient than their 
Bayesian counterparts for this class of models. 
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The multiple imputation estimate Omi is typically obtained by approxi- 
mating the integrals in (3) and (4) using a finite sample from the imputation 
distribution. The procedure can be summarized as follows: 

1. Draw K different realizations (imputations) ~ [S|S] from the 
imputation distribution. 

* (k) - (fc) 

2. For each realization, find the MLE 9 and asymptotic variance 

of the estimate based on the the full data: (S). 

3. Combine results from different imputations using finite sample versions 
of the conditional expectation (3) and variance (4) results: 

^ = ^E^ (5) 
k=i 

Var(6 MI ) = 1 f; var(6 {k) ) + 1 £ (*<*> -e)\ (6) 
k=l k=l 

Equations (5) and (6) are the commonly used combining rules (Rubin, 
1987) for multiple imputation estimators in the scalar case. When is 
vector- valued, similar approaches can be used (Rubin, 1987). 

3 Continuous-Time Discrete-Space Movement Model 

Having described the multiple imputation framework, we now focus on spec- 
ifying a model of animal response to drivers of movement [(S, S)|0] that is 
flexible and computationally efficient. In doing so, we will assume a discrete 
(e.g., gridded) model for space (e.g., Hooten et al., 2010), and model animal 
movement as a continuous-time random walk through the discrete, gridded 
space. 

Let the the study area be defined as a graph (G,A) of M locations 
G = (G±, G2, • • • , Gm) connected by "edges" A = {Ay : i ~ j, i = 1, . . . , M} 
where i ~ j means that the nodes G{ and Gj are directly connected. For 
example, in a gridded space each grid cell is a node and the edges connect 
each grid cell to its first-order neighbors (e.g., a rook's neighborhood). In 
typical studies, the spatial resolution of the grid cells in G will be determined 
by the resolution at which environmental covariates that may drive animal 
movement and selection are available. 
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Figure 1: Continuous-time continuous-space and continuous-time discrete- 
space representations of an animal's movement path. 




A path realization S from the CTCRW model is continuous in time 
and space (Figure 1). If we consider a discrete, gridded space G, then the 
continuous-time, continuous-space path S is represented by a continuous- 
time, discrete-space path (g, r) consisting of a sequence of grid cells g = 
(G^, Gi 2 , . . . , Gi T ) transversed by the animal's continuous-space path and 
the residence times r = (t±, T2, . . • , tt) in each grid cell. 

In practice, this transformation from continuous space to discrete space 
results in a compression of the data to a temporal scale that is relevant to the 
resolution of the environmental covariates that may be driving movement 
and selection. For example, if an animal is moving slow relative to the time 
it takes to traverse a grid cell in G, then the quasi-continuous path S may 
contain a long sequence of locations within one grid cell. The discrete-space 
representation of the path represents this long sequence of locations as one 
observation (a grid cell Ga and residence time Tt). This data compression 
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is especially relevant for telemetry data, in which observation windows can 
span years or even decades for some animals. 



3.1 Random Walk Model 

The discrete-space representation (g, r) of the movement path allows us to 
use standard discrete-space random walk models to make inference about 
possible drivers of movement. While we will relax this assumption later to 
account for temporal autocorrelation in movement behavior, we initially as- 
sume that the the i-th observation (Gi t ,Tt) in the sequence is independent 
of all other observations in the sequence. Under this assumption, the likeli- 
hood of the sequence of transitions {(Gi t — > Gi t+1 , Tt), t = 1, 2, . . . , T} is just 
the product of the likelihoods of each individual observation. We will focus 
on modeling each transition (Gi t — > Gi t+1 ,r t ), and will drop the t subscript 
in this section to simplify notation. 

If an animal is in cell Gi at time t, then define the Poisson rate of 
transition from cell Gi to a neighboring cell Gj as 

Xiji/3) = exp{x^./3} (7) 

where Xjj is a vector containing covariates related to drivers of movement 
specific to cells i and j, and /3 is a vector of parameters that define how 
each of the covariates in are correlated with animal movement. The 
total transition rate Aj from cell i is the sum of the transition rates to all 
neighboring cells: Aj(/3) = Ylj^i^ijiP 1 ) an d the time r that the animal 
resides in cell Gi is exponentially-distributed with rate parameter equal to 
the total transition rate Xi(f3): 

[r| / 3]=A l ( / 3)exp{-rA l (/3)}. (8) 

When the animal transitions from cell Gi to one of its neighbors, the 
probability of transitioning to cell Gk, an event we denote as Gi — > Gk, 
follows a multinomial distribution with probability proportional to the tran- 
sition rate Xik to cell 

If, as is commonly assumed in the study of Markov random walks, the 
residence time and eventual destination are independent events, then the 
likelihood of the observation (Gi — > Gk,r) is the product of the likelihoods 
of its parts: 
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[d -> G k ,r\P] = • A,(/3)exp{-rA,(/3)} 

= A lfc (/3)exp{-rA,(/3)}. (10) 

3.2 Latent Variable Representation 

We now introduce a latent variable representation of the transition process 
that is equivalent to (10), but allows for inference within a standard gen- 
eralized linear modeling framework. For each j such that i ~ j, define Zj 
as 

_ Jl , Gi^Gj 
, o.w. 



and let 



[^rl^ocAgexpl-rA^)} (11) 

where rii is the number of neighbors of the z-th grid cell. Then the product 
of [zj, t\(3\ over all j such that i ~ j is proportional to the likelihood (10) of 
the observed transition: 



n^rlflcx Agexp{-rA y (/3)} 

= A ifc (/3)exp{-rAi(/3)} , where G { G k 
= [G l ^G k ,r\/3] 

The benefit of this latent variable representation is that the likelihood 
of Zj,r\f3 in (11) is equivalent to the kernel of the likelihood in a Poisson 
regression with the canonical log link, where Zj are the observations and 
log(r) is an offset or exposure term. The likelihood of the entire continuous- 
time, discrete-space path (g, r) can be written (returning t to the notation) 
as: 

T 

[g,r|/3] = [Z,r|/3] oc n II [Ajj t (/3)exp{-r t A itjt (/3)}] (12) 

where Z = (zi,...,zy)' is a vector containing the latent variables Zj = 
(zi i: zi 2 ,. . . , Zi K )' for each grid cell in the discrete-space path. Inference can 
be made on (3 in (12) using standard Poisson GLM approaches (e.g., max- 
imum likelihood). This provides a computationally efficient framework for 
the statistical analysis of potential drivers of movement within the multiple 
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imputation framework of Section 2.2. Multiple path realizations (imputa- 
tions) can be drawn from [S|S]. Each continuous path S can then be trans- 
formed into a CTDS path (g, r), which can then be used to make inference 
on (3 using (12). The results from the multiple imputed paths can then be 
combined using (5) and (6), resulting in a multiple imputation estimate (3 MI 
and estimated variance Var((3 MI ). 

4 Time- Varying Behavior and Variable Selection 

In this section we describe how covariate effects can be allowed to vary over 
time using a varying-coefficient model, and how variable selection can be 
accomplished through shrinkage estimation. 

4.1 Changing Behavior Over Time 

Animal behavior and response to drivers of movement can change signif- 
icantly over time. These changes can be driven by external factors such 
as changing seasons (e.g., Grovenburg et al., 2009) or predator/prey inter- 
actions (e.g., Lima, 2002), or by internal factors such as internal energy 
levels (e.g., Nathan et al., 2008). The most common approach to modeling 
time- varying behavior in animal movement is state-space modeling, typically 
within a Bayesian framework (e.g., Morales et al., 2004; Jonsen et al., 2005; 
Getz and Saltz, 2008; Nathan et al., 2008; Forester et al., 2009; Gurarie 
et al., 2009; Merrill et al., 2010). Often, the animal is assumed to exhibit 
a number of behavioral states, each characterized by a distinct pattern of 
movement or response to drivers of movement. The number of states can 
be either known and specified in advance by the researcher (e.g., Morales 
et al., 2004; Jonsen et al., 2005) or allowed to be random (e.g., Hanks et al., 
2011; McClintock et al., 2012). 

State-space models are an intuitive approach to modeling changing be- 
havior over time, but there are limits to the complexity that can be mod- 
eled using this approach. Allowing the number of states to be unknown 
and random requires a Bayesian approach with a changing parameter space. 
This is typically implemented using reversible- jump MCMC methods (e.g., 
Green, 1995; McClintock et al., 2012; Hanks et al., 2011), which are compu- 
tationally expensive and can be difficult to tune. Our approach is to use a 
computationally efficient GLM (12) to analyze parameters related to drivers 
of animal movement. Instead of using the common state-space approach, 
we employ varying-coefficient models (e.g., Hastie and Tibshirani, 1993) to 
model time- varying behavior in animal movement. 
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For simplicity in notation, consider the case where there is only one 
covariate x in the model (7) and no intercept term. The model for the Pois- 
son transition rate (7) will typically contain an intercept term and multiple 
covariates {x}, and the varying-coefficient approach we present generalizes 
easily to this case. In a time-varying-coefficient model, we allow the pa- 
rameter f3(t) to vary over time in a functional (continuous) fashion. The 
transition rate (7) then becomes: 

\ ij {/3(t)) = ex P {x ij /3(t)} 

where t is the time of the observation and Xij is the value of the covariate 
related to the Poisson rate of moving from cell i to cell j. The functional 
regressor (3(t) is typically modeled as a linear combination of n sp i spline 
basis functions {<p k (t), k = 1, . . . , n sp /} : 

^) = £a^(t). 
k=i 

B-spline basis functions are among the most-widely used choices for {4>k(t)}, 
and are appropriate in most cases. Fourier basis functions are commonly 
used for {^fc(i)} when the behavior is thought to be periodic (e.g., diurnal 
variation in behavior). 

Under this varying-coefficient specification, (7) can be rewritten as 

\ ij (P(t)) = exp{xi j 0(t)} 

= exp{<^.a} (13) 

where a = (ai, . . . , a nspl )' and ^ = x^ ■ (0i(i), . . . , 4>n spl (*))'■ The result is 
that the varying-coefficient model can be written as a GLM with a modified 
design matrix. This specification provides a flexible framework for allowing 
the effect of a driver of movement (x) to vary over time that is computa- 
tionally efficient and simple to implement using standard GLM software. 

4.2 Shrinkage Estimation 

The model we have specified in (12) is likely to be overparameterized, es- 
pecially if we utilize a varying-coefficient model (13). Animal movement 
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behavior is complex, and a typical study could entail a large number of 
potential drivers of movement, but an animal's response to each of those 
drivers of movement is likely to change over time, with only a few drivers 
being relevant at any one time. Under these assumptions, many of the pa- 
rameters ak in (13) are likely to be very small or zero. Multicollinearity is 
also a potential problem, as many potential drivers of movement could be 
correlated with each other. 

We propose a shrinkage estimator of a using a lasso penalty (Tibshirani, 
1996). The typical maximum likelihood estimate of a. is obtained by maxi- 
mizing the likelihood [Z,t|o:] from (12), or equivalently by maximizing the 
log-likelihood log[Z, r\a]. The lasso estimate is obtained by maximizing the 
penalized log-likelihood, where the penalty is the sum of the absolute values 
of the regression parameters 




As the tuning parameter 7 increases, the absolute values of the regression 
parameters are "shrunk" to zero, with the parameters that best de- 

scribe the variation in the data being shrunk more slowly than parameters 
that do not. Cross-validation is typically used to set the tuning parameter 
7 at a level that optimizes the model's predictive power. 

Shrinkage approaches such as the lasso are well-developed for GLMs, and 
computationally-efficient methods are available for fitting GLMs to data 
(e.g., Friedman et al., 2010). Recent work has also applied the lasso to 
multiple imputation estimators (e.g., Chen and Wang, 2011). The main 
challenge in applying the lasso to multiple imputation is that a parameter 
may be shrunk to zero in the analysis of one imputation but not in the 
analysis of another. The solution is to use a so-called group lasso (Yuan and 
Lin, 2006), in which a group of parameters is constrained to either all equal 
zero or all be non-zero together. 

The lasso can be seen as a constrained optimization problem, with 
ai asso minimizing the mean squared error subject to the constraint that 
||c*iasso||i < where || • || is the L-l norm and the value of v is determined 
by the tuning parameter 7. As the estimate aq a sso typically falls on the 
boundary of the constrained space, conventional approaches for estimating 
the standard error of ai asso are unavailable. Bayesian versions of the lasso 
(Park and Casella, 2008) and group lasso (Raman et al., 2009) provide al- 
ternatives that allow for estimation of the uncertainty about the parameters 
aiasso through posterior analysis. 
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The Bayesian approach entails significantly more computational com- 
plexity, and may not be as accessible to practitioners. We instead focus on 
the likelihood-based stacked lasso estimate of Chen and Wang (2011). In this 
estimate, instead of computing the lasso estimate ctq asso for each imputation 
individually, and then combining the results using (5) and (6), the imputed 
data from all estimates are "stacked" together and a group lasso estimate is 
obtained for the combined data. We note that this likelihood-based stacked 
lasso approach does not allow for the estimation of the variance of ai asso . If 
uncertainty estimates are a priority, we recommend choosing a parsimonious 
selection of potential drivers of movement a priori that exhibit little mul- 
ticollinearity and computing the traditional multiple imputation estimates 

OLMI- 

5 Drivers of Animal Movement 

We now provide some examples showing how a range of hypothesized drivers 
of movement could be modeled within the CTDS framework. Following 
Hooten et al. (2010), we consider two distinct categories for drivers of move- 
ment from cell Gi to cell Gf location-based drivers ({pki, k = 1, 2, ... , K}) 
which are determined only by the characteristics of cell Gi, and directional 
drivers ({quj, I = 1,2, ... , L}) which vary with direction of movement. Un- 
der a time- varying coefficient model for each driver, the transition rate (7) 
from cell Gi to cell Gj is 

A,, (P(t)) = exp lpo(t) + J2PkA(t) +X>«#(*)| (15) 
I k=i i=i J 

where /3o(t) is a time- varying intercept term, {/?&(£)} are time- varying ef- 
fects related to location-based drivers of movement, and {Pi(t)} are time- 
varying effects related to directional drivers of movement. We consider both 
location-based and directional drivers in what follows. 

5.1 Location-Based Drivers of Movement 

Hooten et al. (2010) denote static, non-directional drivers of movement as 
location-based drivers of movement. Location-Based drivers of movement 
can be used to examine differences in animal movement rates that can be 
explained by the environment an animal resides in. In the CTDS context, 
location-based drivers would be covariates dependent only on the charac- 
teristics of the cell where the animal is currently located. Large positive 
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(negative) values of the corresponding (3k(t) would indicate that the animal 
tends to transition quickly (slowly) from a cell containing the cover type in 
question. 

5.2 Directional Drivers of Movement 

In contrast to location-based drivers, which describe the effect that the local 
environment in which the animal resides has on movement rates, directional 
drivers of movement (Brillinger et al., 2001; Hooten et al., 2010; Hanks et al., 
2011) capture directional selection by the individual. 

A directional driver of movement is defined by a vector which points 
toward (or away) from something that is hypothesized to attract (or repel) 
the animal in question. Let V; be the vector corresponding to the l-th 
directional driver of movement. In the CTDS model for animal movement, 
the animal can only transition from cell Gi to one of its neighbors Gj : j ~ i. 
Let Wjj be a unit vector pointing from the center of cell Gi in the direction 
of the center of cell Gj. Then the covariate quj relating the Z-th directional 
driver of movement to the transition rate from cell Gi to cell Gj is the dot 
(or inner) product of v; and Wj.,-: 

q Uj = vjw^. 

Then puj will be positive when v; points nearly in the direction of cell 
Gj, negative when v; points directly away from cell Gj, and zero if is 
perpendicular to the direction from cell Gi to cell Gj. 

5.3 Examples 

We now provide multiple examples of drivers of movement to illustrate the 
range of effects that can be modeled using this framework. 

5.3.1 Overall Movement Rate 

The intercept term fio{t) in (15) can be seen as a driver of movement in which 
Poi = 1 for every cell G; £ G. This intercept term controls the animal's 
overall rate of transition from any cell, and thus models the animal's overall 
movement rate. Allowing the intercept parameter (3o(t) to vary over time 
could reveal changes in activity levels over time. For example, we might 
expect (3o(t) to be larger at night for nocturnal species and smaller during 
the day. 
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5.3.2 Movement Response to Land Cover Types 

Indicator variables could be used to examine how animal movement differs 
between different landscape cover types (e.g., forest vs. plains) by setting 
Pki = 1 for each cell G% that is classified as containing the fc-th cover type. As 
in the case of the intercept, allowing the parameter (3k(t) related to the k-th 
cover type to vary over time can reveal variation in an animal's movement 
pattern through the cover type. For example, an animal may move quickly 
through open terrain during the day, but may move more slowly through 
the same terrain at night. 

5.3.3 Environmental Gradients 

Animals may use environmental gradients for navigation. For example, a 
mule deer might move predominantly in the direction of increasing elevation 
during a spring migration (e.g., Hooten et al., 2010), or a seal might follow 
gradients in sea surface temperature to navigate toward land (e.g., Hanks 
et al., 2011). Such effects can be modeled by including a directional driver 
of movement in (15) defined by a gradient vector v; which points from the 
center of cell Gi in the direction of steepest increase in the covariate x\. 
Positive values of fii indicate that the animal moves generally towards cells 
with higher values of x\ , while negative values of Pi indicate that the animal 
moves generally towards cells with lower values of x\. 

5.3.4 Activity Centers 

Many animals exhibit movement patterns that are centered on a location 
in space. This central location may be temporary, such as a kill site for a 
predator (e.g., Knopff et al., 2009), or more permanent, such as a den for 
a central place forager (e.g., Hanks et al., 2011; McClintock et al., 2012). 
The relatively new class of spatial capture-recapture models (e.g., Royle and 
Young, 2008) model detection probability as a decreasing function of dis- 
tance from a central location (e.g., the "center" of an animal's home range). 
Movement around an activity center can be modeled in the CTDS framework 
by including a directional driver of movement in (15) defined by a vector v; 
which points from the center of cell Gi to the location of the activity cen- 
ter. Then a positive value for would indicate that the animal is generally 
drawn toward this activity center. If the activity center is considered to be 
temporary (such as a kill site for a predator), then a time-varying-coefficient 
model should be used. The variable selection obtained through the lasso es- 
timate can indicate the range of time in which the animal's movement is 
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centered around the activity center. If the activity center is considered to 
be permanent through the duration of the study, a varying- coefficient model 
may not be needed. 

Under the likelihood-based specification of the CTDS model for animal 
movement, it is necessary to specify the locations of all hypothesized activity 
centers beforehand. In Section 6, we show an example of the specification of 
hypothesized activity centers (potential kill sites for mountain lions) using 
the original telemetry data. If a Bayesian formulation of the CTDS model 
were used, then the location of hypothesized activity centers could be ran- 
dom, and inference could be made on their locations jointly with inference 
on the movement parameters, similar to what is done in spatial capture- 
recapture models (e.g., Royle and Young, 2008). 

5.3.5 Conspecific Interaction 

An animal's movement patterns can be greatly affected by interaction with 
conspecifics. For example, one animal could follow the trail left by another 
animal, two animals could avoid one another by changing course when they 
become close enough to sense the other animal, or a pair of animals could 
maintain proximity as they move together across the landscape. While there 
are many possible approaches to modeling such dependence in behavior, we 
choose to model each of these interactions through the inclusion of direc- 
tional effects in the CTDS modeling framework. For example, a directional 
driver could be included in the movement model for one animal that is 
defined by a vector pointing to the current location of another animal to ex- 
amine whether the animal being modeled is attracted to (/3; > 0) or avoids 
(Pi < 0) the conspecific. 

5.3.6 Directional Persistence 

The CTCRW model of Johnson et al. (2008) is based on a correlated random 
walk model for velocity that allows for directional persistence in animal 
movement. So far, we have assumed that each discrete movement step in 
our CTDS model is independent, but this assumption is not met if the animal 
exhibits any directional persistence. To account for directional persistence 
in the CTDS approach, we use an autoregressive approach by including a 
directional driver of movement at each discrete movement step that is defined 
by a vector pointing in the direction of the previous move. For example, if 
the animal moved west in the previous discrete movement step, then the 
autoregressive vector for the next step points west as well. Positive values 
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of the fi related to this directional driver of movement indicate that the 
animal is likely to maintain its direction of movement over time. 

5.4 Spatial and Temporal Scale 

The choice of scale for a study can greatly influence results (e.g., Boyce, 
2006). When speaking of the scale of a study, one could look at the grain, 
or resolution, at which the process is modeled, or the extent (coverage) over 
which the process is modeled. The spatial and temporal extent of a study of 
animal movement are determined by the telemetry data and the posterior 
predictive path distribution [S|S]. However, when implementing the CTDS 
approach, the researcher must make three choices pertaining to the grain 
or resolution: (1) the temporal scale at which the CTCRW movement path 
of the animal is sampled, (2) the spatial scale of the grid over which the 
discrete-space movement will be modeled, and (3) the temporal scale of the 
varying coefficient model, which is determined by the number and resolution 
of spline knots in the spline basis expansion. 

As the CTCRW model of Johnson et al. (2008) is a continuous-time 
model, we recommend sampling from the movement path at as fine an inter- 
val as is feasible. In practice, this will be limited by computational resources 
and the size of the study. The temporal resolution needs to be fine enough 
that realizations from the posterior predictive path distribution [S|S] are 
quasi-continuous and adequately capture the residence time r in each grid 
cell in the CTDS representation of the movement path. 

The choice of spatial resolution of the raster grid on which the CTDS 
process occurs implicitly specifies a time scale at which the movement pro- 
cess is modeled. Coarser spatial resolution (larger grid cells) will correspond 
to longer residence times r in the CTDS mode. The spatial resolution should 
be chosen so that the time scale at which an animal transitions from one 
grid cell to another is a time scale at which the animal in question can make 
meaningful choices about movement and resource selection. The time scale 
implicit in the choice of spatial resolution can be examined by plotting a his- 
togram of the residence times in the CTDS representation of the movement 
path. 

If the lasso penalty is used, then it is common to choose a saturated 
spline basis expansion in the varying-coefficient model, where one spline 
knot is specified at each data point in time. We recommend specifying 
a spline basis expansion with knots at a similar temporal resolution to the 
temporal resolution of the telemetry data. The lasso penalization will shrink 
the overparameterized expansion to a more parsimonious model that best 
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fits the data. While a finer temporal resolution could be used, the posterior 
predictive path distribution is unlikely to show changes in behavior at time 
scales smaller than the time scale of the original data. Using a coarser 
temporal resolution will force /3(t) to be smooth. This would imply that 
changes in animal behavior are gradual and occur at time scales larger than 
the time scale of the data. 

6 Example: Mountain Lions in Colorado 

We illustrate our CTDS random walk approach to modeling animal move- 
ment through a study of mountain lions (Puma concolor) in Colorado, USA. 
As part of a larger study, a female mountain lion, designated AF79, and her 
subadult cub, designated AM80, were fitted with global positioning system 
(GPS) collars set to transmit location data every 3 hours. We analyze the 
location data S from one week of location information for these two animals 
(Figure 2). 

We fit the CTCRW model of Johnson et al. (2008) to both animals' loca- 
tion data using the 'crawl' package (Johnson, 2011) in the R statistical com- 
puting environment (R Development Core Team, 2012). Ten imputations 
from the posterior distribution of the quasi-continuous path distribution 
[S|S] were obtained at one-minute intervals. The result is a quasi-continuous 
path at extremely fine temporal resolution for each imputation. 

For covariate data, we used a landcover map of the state of Colorado cre- 
ated by the Colorado Vegetation Classification Project (http:/ /ndis. nrel.colostate.edu/coveg/), 
which is a joint project of the Bureau of Land Management and the Col- 
orado Division of Wildlife. The landcover map contained gridded landcover 
information at 100m square resolution. Figure 3 shows a histogram of the 
residence times r in each grid cell in the CTDS representation of the move- 
ment path of AM80. The area traveled by the two animals in our study 
was predominantly forested. To assess how the animals' movement differed 
when in terrain other than forest, we created an indicator covariate where all 
forested grid cells were assigned a value of zero, and all cells containing other 
cover types, including developed land, bare ground, grassland, and shrubby 
terrain, were assigned a value of one (Figure 2a). This covariate and an in- 
tercept covariate were used as location-based covariates in the CTDS model 
for both animals. 

For each animal, we created a set of potential kill sites (PKS) by examin- 
ing the original GPS location data (Figure 2). Knopff et al. (2009) classified 
a location as a PKS if two or more GPS locations were found within 200m 
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Figure 2: Telemetry data for a female mountain lion (AF79) and her male 
cub (AM80). A location-based covariate was defined by landcover that was 
not predominanty forested (a). Potential kill sites were identified, and a 
directional covariate defined by a vector pointing toward the closest kill site 
(b) was also used in the CTDS model. 

(a) Static Cover - Not Forest 









4- Telemetry Locations -AF79 
o Telemetry Locations -AM80 

Forested Terrain 
■ Other Terrain 


















4-*'-'' 







456000 458000 460000 462000 



(b) Distance to Nearest Potential Kill Site For AM80 
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Figure 3: Histogram of residence times in each lOOm-square grid cell in 
the continuous-time discrete-space representation of the movement path of 
a male mountain lion (AM80). 

Histogram of Transition Times (t) for AM80 
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of the site within a six-day period. We added an additional restraint that 
at least one of the GPS locations be during nighttime hours (9pm to 6am) 
for the point to be classified a PKS. We then created a covariate raster layer 
containing the distance to the nearest PKS for each grid cell (Figure 2). A 
directional covariate defined by a vector pointing towards the nearest PKS 
was included in the CTDS model for both animals. 

To examine how the movement path of the cub AM80 affected the move- 
ment path of the mother AF79, we included a directional covariate in the 
CTDS model for AF79 defined by a vector pointing from the mother's loca- 
tion to the cub's location at each time point. Similarly, we included a direc- 
tional covariate in the CTDS model for AM80 defined by a vector pointing 
from the cub's location to the mother's location at each time point. 

For each animal, we also included a directional covariate pointing in the 
direction of the most recent movement at each time point. This covariate 
measures the strength of correlation between moves, and thus the strength of 
the directional persistence shown by the animal's discrete-space movement 
path. As we are assuming an underlying correlated movement model (the 
CTCRW model of Johnson et al. (2008)), we expect the CTDS movement 
to be correlated in time as well. 

To allow for varying behavior over time, we used a varying-coefficient 
model for each covariate in the model. For all covariates except the direc- 
tional covariate related to directional persistence, we used a B-spline basis 
expansion, with regularly spaced spline knots at 3-hour intervals. 

We fit the CTDS model for each path using the 'glmnet' R package 
(Friedman et al., 2010), using a Lasso penalty, with tuning parameter chosen 
to minimize the average squared error of the fit in a 10-fold cross-validation. 

6.1 Results 

The time-varying results for the location-based and directional drivers of 
movement for AF79 are shown in Figure 4, with the corresponding results 
for AM80 shown in Figure 5. As we used a lasso penalty, we can only obtain 
point estimates (confidence intervals are unavailable) of the time-varying 
effects {/3(t)}. A comparison of the differences between the results for AF79 
and AM80 yields some insight into how the movement patterns of these two 
animals differ. 

The intercept effect measures the animal's general movement rate over 
time. Figure 4(a) and Figure 5(a) show the time- varying deviation from the 
grand mean for each animal. For example, the male AM80 moved at a higher 
than average rate of speed near sunrise on Julian day 90, as evidenced by the 
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Figure 4: Time- varying results for the location-based and directional covari- 
ates in the continuous-time discrete-space model for a female mountain lion 
(AF79). Shaded regions indicate nighttime hours. 
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Figure 5: Time- varying results for the location-based and directional covari- 
ates in the continuous-time discrete-space model for a male mountain lion 
(AM80). Shaded regions indicate nighttime hours. 

(a) Static Intercept (Deviation tram Overall Intercept = 3.66) (b) Static Cover - Not Forest 
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positive P(t) in Figure 5(a). The time- varying intercept for both lions tends 
to be higher during nighttime hours and crepuscular periods than during the 
day, indicating higher overall rates of movement during nocturnal hours. 

The location-based response to non-forested land cover (Figure 4b and 
Figure 5b) is zero for much of the study window for both AF79 and AM80. 
Forested land cover makes up over 90% of the study area, thus the moun- 
tain lions encounter different land cover infrequently. On day 90, the male 
(AM80) moves through a patch of non-forested land cover at a relatively 
high rate of speed. On day 89, the female moves through a patch of non- 
forested land cover at a low rate of speed. Both animals encounter non- 
forested patches at other times during the observation window, but the f3(t) 
related to movement through non-forested terrain is near-zero during these 
times. This f3{t) has been shrunk to zero by the Lasso procedure during 
these times, indicating that during these times the animal's movement is 
not greatly different in non-forested patches than it is in forested patches. 

The male cub (AM80) has a fairly consistent positive response to the 
directional covariate pointing toward the nearest PKS from day 88 to day 90 
(Figure 5c), indicating that much of the male's movement during these three 
days resembles a random walk centered on an attractive central location (a 
PKS). During day 91, the negative response to the direction to the nearest 
PKS indicates that the animal is moving away from the nearest PKS. At 
the same time, the positive response to the directional covariate pointing 
towards AF79 (Figure 5d) indicates that AM80 is moving towards AF79 
during this time period. 

Positive values of the autoregressive parameter (Figure 4e and Figure 5e) 
indicate that the animals are in a state of correlated movement (directional 
persistence). The magnitude of the autoregressive parameter is greater in 
general for the female (Figure 4e) than for the male (Figure 5e), indicat- 
ing that the female may have a greater tendency for directed movements 
than the male. The CTDS movement paths are based on the CTCRW 
model of animal movement, which results in correlated movement paths. 
The inclusion of this autoregressive parameter is meant to account for the 
correlation in the CTDS path representation of the underlying correlated 
CTCRW movement path. 

7 Discussion 

While we have couched our CTDS approach in terms of modeling animal 
movement, we can also view this approach in terms of resource selection 
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(e.g., Manly, 2002). Johnson et al. (2008) describe a general framework 
for the analysis of resource selection from telemetry data using a weighted 
distribution approach where an observed distribution of resource use is seen 
as a re-weighted version of a distribution of available resources, and the 
resource selection function (RSF) defines the preferential use of resources by 
the animal. Warton and Shepherd (2010) and Aarts et al. (2012) describe a 
point-process approach to resource selection that can be fit using a Poisson 
GLM, similar to the CTDS model we describe here. In the context of Warton 
and Shepherd (2010), the CTDS approach can be viewed as a resource- 
selection analysis with the available resources at any time defined as the 
neighboring grid cells. The transition rate (15) of the CTDS process to each 
neighboring cell contains information about the availability of each cell, as 
well as the RSF that defines preferential use of the resources in each cell. 

It is notable that the entire analysis in Section 6 requires less than ten 
minutes using a computer with 4 GB of memory and a 1.67 GHz quad- 
core processor. This increase in computational efficiency relative to the ap- 
proaches of Johnson et al. (2008), Hooten et al. (2010), Hanks et al. (2011), 
and McClintock et al. (2012) allows for inference on complex behavior at 
finer temporal resolution than has been possible previously. To make our 
CTDS approach accessible to practitioners, we have created an R-package 
('ctds') that contains R-code to fit the CTDS model using multiple imputa- 
tion as described in Sections 2-4. A script file contained in the 'ctds' package 
allows for the re-analysis of the telemetry data of the two mountain lions 
analyzed in Section 6. This R-package can be downloaded from the first 
author's website (http://www.stat.colostate.edu/~hanks). 

The CTDS approach to modeling animal movement is flexible, and can 
be extended using standard approaches for GLMs. For example, if population- 
level inference is desired, the movement paths from multiple animals could be 
analyzed jointly, with population-level parameters in the GLM being shared 
by all animals. Similarly, interaction terms could be included in the model 
by including multiplied covariates in the design matrix. Fitting movement 
models in a GLM framework allows for many natural extensions with little 
additional effort. 

We have focused on building an individual model for animal movement 
that is intuitive and computationally-efficient. If population- level inference 
on data from multiple animals is desired, there are multiple potential ap- 
proaches. The first is to analyze movement paths from multiple animals 
jointly, with population-level parameters in the GLM being shared by all 
animals and individual variation modeled using standard random effects ap- 
proaches. However, this approach may not be straightforward to implement 
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or interpret in the context of a varying-coefficient model. Each individual 
animal is likely to encounter and be influenced by different drivers of move- 
ment (e.g., local environmental factors or nearby conspecifics) at differing 
times throughout the observation window. In this situation, it may make 
little sense to examine a population-level response to a particular driver of 
movement at a particular time, as in a typical random-effects analysis. 

Instead, one approach to a population-level analysis could be a post-hoc 
analysis of the time- varying response to covariates (3(t). For example, Hanks 
et al. (2011) used a cluster analysis to examine different movement regimes 
shared across individuals in the population. Differences in the movement 
patterns exhibited by subgroups (e.g., male vs. female) were then examined 
by a comparison of the proportion of time spent by each subgroup in each 
of the movement regimes. 

The use of directional drivers of movement has a long history. Brillinger 
et al. (2001) model animal movement as a continuous-time, continuous-space 
random walk where the drift term is the gradient of a "potential function" 
that defines an animal's external drivers of movement. Tracey et al. (2005) 
use circular distributions to model how an animal moves in response to a vec- 
tor pointing towards an object that may attract or repel the animal. Hanks 
et al. (2011) and McClintock et al. (2012) make extensive use of gradients 
to model directed movements, and movements about a central location. In 
our study of mountain lion movement data, we used directional drivers of 
movement to model conspecific interaction between a mother (AF79) and 
her cub (AM80). Interactions between predators and prey could also be 
modeled using directional covariates defined by vectors pointing between 
animals. Some movements based on memory could also be modeled using 
directional covariates. For example, a directional covariate defined by a vec- 
tor pointing to the animal's location one year prior could be used to model 
seasonal migratory behavior. The ability to model both location-based and 
directional drivers of movement make the CTDS framework a flexible and 
extensible framework for modeling complex behavior in animal movement. 
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